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Abstract 

Implicit  integration  schemes,  such  as  Runge-Kutta  and  Runge-Kutta-Nystrom  meth¬ 
ods,  are  widely  used  in  mathematics  and  engineering  to  numerically  solve  ordinary 
differential  equations.  Every  integration  method  requires  one  to  choose  a  step-size, 
h,  for  the  integration.  If  h  is  too  large  or  too  small  the  efficiency  of  an  implicit  scheme 
is  relatively  low.  As  every  implicit  integration  scheme  has  a  global  error  inherent 
to  the  scheme,  we  choose  the  total  number  of  computations  in  order  to  achieve 
a  prescribed  global  error  as  a  measure  of  efficiency  of  the  integration  scheme.  In 
this  paper,  we  propose  the  idea  of  choosing  h  by  minimizing  an  efficiency  function 
for  general  Runge-Kutta  and  Runge-Kutta-Nystrom  integration  routines.  This  effi¬ 
ciency  function  is  the  critical  component  in  making  these  methods  variable  step-size 
methods.  We  also  investigate  solving  the  intermediate  stage  values  of  these  routines 
using  both  Newton’s  method  and  Picard  iteration.  We  then  show  the  efficacy  of  this 
approach  on  some  standard  problems  found  in  the  literature. 

Key  words:  Runge-Kutta,  Runge-Kutta-Nystrom,  numerical  integration,  variable 
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1  Introduction 

Recently,  there  has  been  interest  in  the  literature  concerning  the  use  of  geo¬ 
metric  integration  methods,  which  are  numerical  methods  that  preserve  some 
geometric  quantities.  For  example,  the  symplectic  area  of  a  Hamiltonian  sys¬ 
tem  is  one  such  concern  in  recent  literature  [1-4] .  Tan  [5]  explores  this  concept 
using  implicit  Runge-Kutta  integrators.  Hamiltonian  systems  are  of  particular 
interest  in  applied  mathematics,  and  in  fact  we  test  our  variable  step-size  selec¬ 
tion  method  on  a  well-known  Hamiltonian  system  in  Section  4.2.  Furthermore, 
Hairer  and  Wanner  [6,7]  showed  that  although  implicit  Runge-Kutta  methods 
can  be  difficult  to  implement,  they  possess  the  strongest  stability  properties. 
These  properties  include  A-stability  and  A-contractivity  (algebraic  stability). 
These  are  the  main  reasons  we  choose  to  investigate  variable  integration  step- 
size  selection  using  Runge-Kutta  methods. 

First  order  ordinary  differential  equations  are  solved  numerically  using  many 
different  integration  routines.  Among  the  most  popular  methods  are  Runge- 
Kutta  methods,  multistep  methods  and  extrapolation  methods.  Hull,  Enright, 
Fellen  and  Sedgwick  [8]  have  written  an  excellent  comparison  of  these  types 
of  methods.  They  test  a  number  of  Runge-Kutta  methods  against  multistep 
methods  based  on  Adams  formulas  and  an  extrapolation  method  due  to  Bu- 
lirsch  and  Stoer  [9] .  A  goal  of  that  paper  was  to  compare  these  different  types 

Email  addresses:  raymond.w.holsapple@ttu.edu  (Raymond  Holsapple), 
ram.iyer@ttu.edu  (Ram  Iyer),  david.doman@wpafb.af.mil  (David  Doman). 
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of  methods  as  to  how  they  handle  routine  integration  steps  under  a  variety 
of  accuracy  requirements.  Implicit  or  explicit  integration  methods  require  one 
to  choose  a  step-size,  h,  for  the  integration.  One  of  the  questions  Bulirsch 
and  Stoer  investigate  is  a  strategy  for  deciding  what  step-size  h  to  use  as  the 
methods  progress  from  one  step  to  another.  Others  have  investigated  this  very 
same  problem  in  the  past  [8,10-12]. 

In  this  paper,  we  propose  the  idea  of  choosing  variable  step-sizes  by  minimiz¬ 
ing  an  efficiency  function  for  general  Runge-Kutta  and  Runge-Kutta-Nystrom 
integration  routines.  As  every  implicit  integration  scheme  has  a  global  error 
inherent  to  the  scheme,  we  choose  the  total  number  of  computations  in  order 
to  achieve  a  prescribed  global  error  as  a  measure  of  efficiency  of  the  integra¬ 
tion  scheme.  For  illustration  purposes,  consider  Figure  1.0.1,  referring  to  the 
solution  of  (2).  Let  x(tk)  be  our  approximation  to  x(tk).  We  determine  the 
variable  step-sizes  hi,  h2, . . . ,  h$,  where  =  4  —  4-i,  so  that  we  minimize  an 
efficiency  function  that  minimizes  the  sum  of  the  total  number  of  computa¬ 
tions  to  compute  x(tk)  for  k  —  1,  2, . . . ,  8  and  the  global  error  that  propagates 
from  the  local  truncation  errors  at  each  step  of  integration.  An  alternate  and 
perhaps  simple  way  of  understanding  our  method  is  that  we  choose  hL||A||, 
where  L  is  a  Lipschitz  constant  for  the  differential  equation  and  A  is  a  matrix 
that  describes  the  integration  scheme  used  to  integrate  the  differential  equa¬ 
tion.  Of  course,  in  a  fixed  step-size  method  one  chooses  h.  To  the  best  of  our 
knowledge,  our  proposed  method  is  novel. 

In  the  rest  of  this  section,  we  briefly  describe  the  approaches  found  in  the 
literature.  Hull,  Enright,  Fellen  and  Sedgwick  [8]  approach  this  topic  as  follows. 
First,  they  determine  hmax ,  which  is  a  measure  of  the  “scale”  of  a  problem.  This 
helps  to  allow  them  from  not  stepping  past  any  interesting  fluctuations  in  the 
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Fig.  1.0.1.  Illustration  of  variable  step-sizes  and  error  propagation  in  numerical 
integration 

solution.  Then,  for  their  Runge-Kutta  methods  they  compute  r,  an  estimate 
on  the  local  truncation  error,  which  must  be  bounded  by  the  tolerance,  s. 
They  then  compute 


hnew  =  min  {hmax,  0.9 holi  (£/t)1/p}  .  (1) 

where  p  is  the  order  of  the  Runge-Kutta  routine  being  used. 

Stoer  and  Bulirsch  [10]  arrive  at  a  very  similar  solution  to  this  problem.  To  de¬ 
scribe  what  they  do,  we  first  note  that  throughout  this  paper,  we  will  consider 
solving  the  following  first  order  ordinary  differential  equation: 

cIt 

—  =  f(t,x),  x(0)  =  x0  £  !Rn,  (2) 

where  /  :Rx  ]Rn  — »•  ]Rn  is  Lipschitz  continuous  in  the  second  argument,  i.e. 
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for  any  f  G  R  and  any  vectors  x  G  IRn,  y  G  !Rn ,  we  have 


II f(t,x)  -  f(t,y) ||  <  L(t)\\x  —  y\\,  (3) 

where  L(-)  G  L°°[0,  T].  Stoer  and  Bulirsch  consider  two  discretization  methods, 
$1  and  d>2,  of  Runge-Kutta  type  to  solve  (2),  the  first  with  order  p  and  the 
second  with  order  p  +  1,  i.e.  we  have 


Xk+i  =  Xk  +  h0]d$i{tk,xk’,hM)  (4) 

xk+1  =  xk  +  hold<F 2 (4 ,  xk ;  hold)  (5) 


Denoting  the  tolerance  by  e,  and  given  a  current  step-size,  hold ,  they  obtain: 


firiRW  hi 


old 


•A;+i  ^fc+i 


i/(p+i) 


(6) 


Stoer  and  Bulirsch  go  on  to  recommend  after  extensive  numerical  experimen¬ 
tation,  that  equation  (6)  be  altered  to 


/i„nw  =  0.9/q 


eh, 


old 


1/p 


(7) 


As  one  can  see,  formulas  (1) ,  (6)  and  (7)  depend  on  some  measure  of  the  local 
error  at  the  (k  +  l)st  step  of  integration.  Stoer  and  Bulirsch  [10]  also  point  out 
that  there  is  another  way  to  determine  hnew,  but  it  requires  one  to  estimate 
higher  order  derivatives  of  f.  For  example,  a  fourth  order  Runge-Kutta  method 
would  require  one  to  estimate  derivatives  of  /  of  the  fourth  order.  Not  only  is 
this  very  costly,  but  this  other  method  uses  the  local  truncation  error  at  the 
fcth  step  of  integration. 


Houwen  [11]  took  a  similar  approach  to  adjusting  step-sizes.  Again  let  £  be 
the  tolerance.  Houwen  then  forms  a  discrepance  function  d(tk,xk ;  hoid)  at  the 
point  ( tk,xk ).  Then  the  new  step-size  is  determined  to  be  the  solution  of  the 
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equation 


\\d(tk,  xk;  hold)\\  =  e.  (8) 

Houwen  considers  three  types  of  discrepance  functions: 

(1)  an  approximation  to  the  local  discretization  error 

(2)  the  residual  term  left  when  the  local  difference  solution  is  submitted  into 
the  differential  equation 

(3)  the  discrepance  of  linearity  of  the  differential  equation 

The  first  two  clearly  are  functions  that  are  some  measure  of  local  error  of  the 
difference  scheme  being  used.  The  discrepance  of  linearity  method  is  merely 
a  way  to  choose  hnew  such  that  the  Jacobian  matrix  for  non-linear  systems 
does  not  change  very  much,(i.e.  within  some  tolerance  s).  over  the  interval 
[tk,tk  +  hnevr] •  This  method  also  deals  with  some  measure  of  local  stability  of 
the  differential  equation. 

Cano  and  Duran  [12]  investigate  variable  step-size  selection  using  linear  mul¬ 
tistep  methods.  Again  consider  equation  (2).  Given  a  tolerance  £,  they  let 

K  =  es(x(tn),e)  +  0(ep),  (9) 

where  p  is  the  order  of  the  method  and  s  is  function  satisfying  the  following: 

(1)  •Smin  ^  s(x ,  £:)  ^  Smax,  with  Smjn,  Smax  >  0, 

(2)  s  is  C°°  in  both  arguments  and  all  the  derivatives  of  s  are  bounded. 
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2  Implicit  Integration  Methods 


Numerical  methods  for  solving  initial  value  problems  such  as  (2)  may  be  either 
explicit  or  implicit.  The  focus  of  this  paper  is  concentrated  on  using  implicit 
methods.  In  this  section,  we  describe  two  classes  of  implicit  numerical  inte¬ 
gration  schemes  and  how  one  might  use  the  methods  to  solve  (2).  We  assume 
the  solution  exists  for  t  £  [0,  T],  with  T  >  0. 


2.1  Runge-Kutta  Methods 

A  detailed  description  of  a  general  s-stage  Runge-Kutta  method  for  the  so¬ 
lution  of  (2)  can  be  found  in  many  publications  [1,3,13].  The  general  method 
for  a  fixed  step-size,  h,  is  described  below: 


yik=xk  +  h^2aijf(tk  +  cjh,yjk),  i  =  l,...,s,  (10) 

3= 1 
s 

xk+i  =  xk  +  h^2bif(tk  +  Cih,yik),  xo  =  x(0).  (11) 

2=1 


In  the  above  equations,  the  yik  are  stage- values  that  must  be  computed  at 
every  step  of  the  integration,  and  Xk  approximates  the  exact  solution  x ( A )  at 
the  point  A  =  kh,  where  h  is  the  fixed  step-size  of  integration.  The  al3  and  bt 
are  unique  to  any  particular  Runge-Kutta  scheme  and  the  <:■,  satisfy 


«  =  i  =  1,  •  •  •  ,S  (12) 

3= 1 
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For  notational  purposes,  define  the  following: 


A  = 


an  <212  • 

^1  s 

yik 

Cl 

Xk 

&21  a22  •  • 

*  &2S 

,  Yk  = 

...  to 

?r 

,  c  = 

C2 

,  Xk  = 

Xk 

&sl  &S2 

&SS 

Vsk 

Cs 

Xk 

A  =  A®  I 


(13) 

where  I  is  the  n  x  n  identity  matrix  and  0  is  the  Kronecker  product.  In 
other  words,  A  is  the  ns  x  ns  matrix  direct  product  of  A  and  I.  Furthermore, 
consider  the  function  /:Rx  IR"5  — ►  fRns  defined  by 

f(tk  +  cih,ylk) 
f(tk  +  c2h,y2k ) 


f(tk,Yk)  = 


(14) 


f(tk  +  csh,ySk) 

Now  we  can  write  the  system  of  ns  equations  given  in  equation  (10)  as 

Yk  =  Xk  +  hAf(tk,Yk).  (15) 


For  each  k  this  is  an  implicit  equation  involving  the  vectors  {yik  }*=1.  Equation 
(15)  can  be  solved  using  Newton’s  method  or  fixed  point  iteration  (Picard 
iteration).  Let’s  consider  Picard  iteration.  We  solve  (15)  for  each  k  using  the 
following  iterative  scheme: 

Yt'  =  Xk  +  hAf  (4,  Y{)  =  F  (4,  Y>)  .  (16) 


For  any  fixed  k,  the  iterative  scheme  given  in  (16)  will  converge  to  the  solution 


of  (15)  provided  that  F  satisfies  a  favorable  condition.  The  following  theorem 
addresses  this  convergence. 

Theorem  2.1  Consider  the  iterative  scheme  given  by  (16).  Let  L(t )  be  the 
function  from  (3),  and  let  A  be  the  s  x  s  matrix  from  (13).  If  hL(tfc)||A||  <  1 
then  there  exists  a  unique  vector  Y  G  Hns  such  that  F(tk,Y)  =  Y  for  any 
point  tk  €  [0,T]  that  is  fixed.  Furthermore,  the  sequence  Yf+l  =  F(tk,Yf) 
converges  linearly  to  Y . 

Proof.  First,  let  us  refer  to  a  very  useful  result  concerning  the  norms  of 
Kronecker  products.  According  to  Van  Loan  [14],  we  get 


A 


=  P®/||  = 


(17) 


Next,  we  show  that  F  is  Lipschitz  continuous  in  the  second  argument  with 
Lipschitz  constant  hL(4)||A||.  We  will  begin  by  finding  the  Lipschitz  constant 
of  the  function  /.  Choose  any  tk  €  [0,  T]  and  any  vectors  u  and  v  in  IRns  where 
u  =  [(«i)T  («2)t  ' ' '  (W)T  and  ui  £  lRn  for  *  =  1, . . . ,  s.  Similarly  form 

the  vector  v. 


f(tk,u )  -  f(tk,v) 


(/(4,«i)  -  f(tk,vi))T  ■ 


<  L(tk) 
=  L(tk) 


(ui  -  ni)T  •  •  •  (us 
\Ul)T  •••  (us)T }T 


=  L(tk ) 


u  —  v 


(. f(tk,us )  -  f(tk,Vs)) 


MT  •••  (^)T 


(18) 


(19) 

(20) 


Now  we  may  compute  the  Lipschitz  constant  of  F.  Again,  choose  any  tk  € 
[0,  T]  and  any  vectors  u,  v  G  MTi's . 
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II F(tk,  u)  -  F(tk,  v )  ||  =  Xk  +  hAf(tk ,  u)-Xk-  hAf(tk ,  v) 


=  h 


A  (f(tk,u)  -  f(tk,v )) 


<  h 


A 


<  hL(tk ) 
=  hL(tk ) 


/(**,«)  -  /(4,v) 


A 

\u  —  v\\ 

IW 

II 

\u  —  v  1 

(21) 

(22) 

(23) 

(24) 

(25) 


This  shows  that  F  is  Lipschitz  continuous  in  the  second  argument  with  Lip- 
schitz  constant  hL(tfc)||A||.  Since  /iA(U,)||A||  <  1,  we  may  apply  the  Contrac¬ 
tive  Mapping  Theorem  to  F,  [15].  Hence,  there  exists  a  unique  point  Y  G  Etns 
such  that  F(tk,Y )  =  Y  for  any  point  tk  G  R  that  is  fixed.  The  theorem  also 
ensures  that  Y  must  be  the  limit  of  every  sequence  obtained  from  (16)  with  a 
starting  point  Ffc°  G  !Rns.  □ 


Theorem  2.1  suggests  how  one  might  implement  equations  (11)  and  (16)  to 
solve  (2)  on  [0,T],  The  starting  vector  :r0  G  Hn  is  known.  In  general,  assume 
xk  is  known.  Use  the  following  procedure  to  compute  xk+1. 


(1)  Choose  a  tolerance  £  >  0  as  small  as  you  wish. 

(2)  Choose  a  starting  guess  for  the  s  stage-values,  and  denote  this  guess  as 

n°- 

(3)  For  j  =  0, 1, 2, ...,  compute  the  following: 


(a)  Y’+' =  F  (tt,Y‘) 

(b)  S  =  ||Vfc+1  -Y£\. 

(4)  If  6  <  e,  let  Yk  =  Y’+1. 

(5)  Use  the  s  n  x  1  stage-value  vectors  determined  in  step  four  to  explicitly 


compute  xk+\. 


The  method  described  above  is  known  as  Picard  iteration;  Newton’s  method 
might  also  be  used  to  solve  for  the  stage- values.  A  theorem  on  the  convergence 
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of  Newton’s  method  is  more  complicated  than  Theorem  2.1;  it  is  not  sufficient 
to  assume  hL(t*,)||^4||  <  1  in  order  to  guarantee  that  Newton’s  method  con¬ 
verges.  The  Newton-Kantorovich  Theorem  [10,16,17]  provides  sufficient  condi¬ 
tions  for  existence  of  a  solution  to  the  iteration  and  the  uniqueness  of  that  so¬ 
lution.  To  solve  for  the  stage- values  using  Newton’s  method,  step  three  should 
be  replaced  by  the  following: 


(3)  Define  G(tk,Yk )  =  F(t^Yk)  —  17,  and  for  j  =  0,1,2,...,  compute  the 
following: 

(a)  Y=+'  =  Yl  -  ( DG  (4,  iff1)  G  (tkX) 

(b)  S  =  \\G(tt,Y’+1) 


where  DG  represents  the  Jacobian  matrix  of  G. 


2.2  Rung e-Kutta-Ny strom  Methods 

Runge-Kutta-Nystrom  (RKN)  methods  are  similar  to  Runge-Kutta  methods, 
but  are  designed  to  solve  second-order  systems.  Consider  the  following  system: 

dp  x 

=  x(0)  =  x0  e  Hn,  i(0)=D06R’1,  (26) 

which  can  be  written  as 

rlv  dr 

-  =  /(t,x);  -  =  v,  x(0)  =  xo,  v(0)  =  v0.  (27) 

We  assume  a  solution  exists  on  [0,  T]  for  T  >  0.  A  general  s-stage  RKN  method 
may  be  used  to  solve  (27)  on  [0,  T],  and  is  described  by  J.M.  Sanz-Serna  and 
M.P.  Calvo  [3]  as  follows: 
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i  =  1,  ...,s 


(28) 


Vik  =  xk  +  h'jiVk  +  h2^2  aijf(tk  +  ijh,  yjk), 

3= 1 


Vk+i  =  Vk  +  h'%2bif(tk+'Yih,yik)  (29) 

2=1 

S 

xk+i  =  xk  +  hvk  +  h 2  Z)A/(4  +  7i^yiJ-  (3°) 

2=1 

Exactly  as  with  the  Runge-Kutta  methods,  the  ylk  are  the  stage-values  and 
must  be  solved  for  implicitly.  In  addition  to  the  definitions  in  (13),  define 


" 

■ 

o 

o 

vk 

o 

to 

,vk  = 

Vk 

:  0 

1 

o 

o 

C/D 

i _ 

Vk 

,r  =  r®  /, 


(31) 


where  I  is  the  n  x  n  identity  matrix.  Now,  the  ns  equations  in  (28)  may  be 
written  as 

Yk  =  Xk  +  hWk  +  h2Af(tk,  Yk).  (32) 


Again,  we  have  an  implicit  equation  for  the  s  stage- value  vectors  {yik} |=1  C 

]Rn. 


Just  as  in  Section  2.1,  we  may  solve  (32)  using  Newton’s  method  or  by  using 
Picard  iteration.  If  we  use  Picard  iteration,  we  have  the  following  iterative 
scheme: 

Yk+l  =Xk  +  hVVk  +  h2Af(tk,  Yi)  =  H(tk,  Yi).  (33) 

Similar  to  the  proof  of  Theorem  2.1,  one  can  easily  prove  the  following  theorem. 

Theorem  2.2  Consider  the  iterative  scheme  given  by  (33).  Let  L(t )  be  the 
function  from  (3),  and  let  A  be  the  s  x  s  matrix  from  (13).  If  h2 L(tk)\\A\\  <  1 
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then  there  exists  a  unique  vector  Y  e  Hns  such  that  H(tk,Y )  =  Y  for  any 
point  tk  €  [0 ,T]  that  is  fixed.  Furthermore,  the  sequence  Yf+1  =  H(tk,Yf) 
converges  linearly  to  Y. 

Proof.  See  the  proof  of  Theorem  2.1.  □ 

If  we  choose  to  use  Newton’s  method  to  solve  (32)  for  the  stage-values,  it 
is  not  sufficient  to  assume  h2L(tk)\\A\\  <  1.  Once  again,  we  may  refer  to 
Stoer  and  Bulirsch  [10]  for  conditions  that  guarantee  convergence  of  Newton’s 
method. 


3  Step-Size  Selection 

When  implementing  a  Runge-Kutta  (or  Runge-Kutta-Nystrom)  numerical  in¬ 
tegration  routine,  we  have  shown  it  is  sufficient  to  assume  that  hL(tfc)||^4||  <  1 
(or  h2L(tk)\\A\\  <  1)  to  guarantee  convergence  of  the  implicit  scheme  when 
using  a  Picard  iteration.  One  might  wonder  though,  is  there  an  optimal  choice, 
in  the  sense  of  computational  efficiency,  for  hi  If  so,  how  might  it  be  found? 

3.1  Optimization  Using  Picard  Iteration 

Consider  solving  (2)  numerically  on  the  interval  [0,  T]  which  is  partitioned  by 
the  following  sequence  of  points:  {kh}^=0.  In  the  k-th  sub-interval  the  conver¬ 
gence  of  the  Picard  iteration  is  linear,  so  the  number  of  computations  required 
for  convergence,  to  within  s,  to  the  fixed  point  of  (16)  can  be  written  as  a  func¬ 
tion  of  the  Lipschitz  constant  of  the  function  F  :  (hL(tk)\\A\\) .  Then 

the  total  number  of  computations  over  the  interval  [0,  T]  can  be  written  as 
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N(h)  =  EfcLi  iVfe.  In  the  following,  we  find  an  explicit  expression  for  </>(•)  for 
Runge-Kutta  methods.  Consider  the  following  inequalities: 


h+1 


Yi 


F(tk,Y’)  -  F^Yt1) 
<hL(tk)\\A\\  | Yi-Yf1 
<(hL(tk)\\A\\)2  Yt'-Y} 


-2 


(34) 

(35) 

(36) 


<{hHtk)\\A\\Y 
=  Ck(hL(tk)\\A\\Y 


U‘  -  Yk° 


(37) 

(38) 


where  Cj.  =  \\Yk  —  Yk \ |  is  fixed  for  each  k  and  depends  on  the  guess  Yk .  Since 


/iL(4)||74||  <  1,  we  must  have 

3=  Yi+1-Yi 


Yi+1  -  Y’ 


0  as  j  — >  oo.  Suppose  we  want 


<  e;  then,  it  is  sufficient  to  have  Ck  ( hL{tk)\\A\\)3  <  e.  As  a 
stopping  criterion  in  the  /c-th  step  of  integration,  we  choose  to  stop  the  fixed 
point  iteration  at  the  smallest  natural  number  greater  than  or  equal  to  M 
where  M  satisfies  Ck  (hL(tk)\\A\\)M  —  e.  Then  we  have 


M  = 
Nk  = 


In  (g/ Ck) 

In  (hL(tk)\\A\\) 
\M] 


(39) 

(40) 


Now  let  C  =  max  Ck  and  L  =  sup  L(t).  In  (39),  £  and  Ck  depend  on  the  user. 

k  te[o,T] 

Once  these  are  chosen,  the  choice  of  h  depends  on  the  differential  equation 
being  solved,  through  the  Lipschitz  constant  L(tk),  and  on  the  integration 
method  being  implemented,  through  ||^4||.  We  will  try  to  minimize  M  by 
adjusting  the  choice  of  h  to  the  problem  parameters  L(tk)  and  ||R|| .  Notice 
that  C(hL\\A\\)M  =  e  implies  that  Ck(hL(tk)\\A\\)M  <  £  for  each  k.  Thus,  we 
minimize  the  cost  function 
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(41) 


In  {e/C) 


Ji(h)  =  K—r±/ 


In  (AiPH) 

Tln{e/C) 

~  hln  (hL\\A\\y 

This  is  the  same  as  maximizing  the  cost  function 

hln  (hL\\A\\) 


Mh)  = 


Tin  {e/C) 


(42) 


(43) 


By  computing  arg  min  J2{h),  one  finds  the  step-size  h  that  minimizes  the  num¬ 
ber  of  computations  for  the  iterative  scheme  to  converge.  If  this  were  the  only 
measure  of  optimality  of  concern,  it  is  easily  shown,  through  a  calculus  argu¬ 
ment,  that  the  cost  function  J2(h)  is  maximized  when 


h  = 


1 

3|RT 


(44) 


However,  one  might  also  want  the  global  error  of  the  numerical  solution  to  be 
as  small  as  possible.  Global  error  in  any  numerical  integration  scheme  depends 
on  the  scheme  being  used.  In  this  paper,  we  are  concentrating  on  Runge-Kutta 
schemes.  The  global  error  for  Runge-Kutta  schemes  also  varies  depending  on 
the  scheme  one  chooses  to  implement.  For  the  purpose  of  explanation,  let  us 
consider  the  implicit  midpoint  rule.  The  implicit  midpoint  rule  is  a  one-stage 
Runge-Kutta  method  where  an  =  \  and  b\  =  1  in  (10)  and  (11).  The  implicit 
midpoint  method  has  global  error  G{Th2).  Then  to  find  h,  we  alter  the  cost 
function  J2  and  maximize  the  following  cost  function: 


Mh) 


/iln(/iL||H||) 

Tln{s/C) 


-  nTh2, 


(45) 


where  n  is  a  number  to  be  chosen.  First  we  compute 


dJ3  _  ln{hL\\A\\)  +  1 
dh  ~  Tln{e/C) 


(46) 
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and 


d2J 3 


1 


dh2  Thln(e/C) 


-  2  kT. 


(47) 


In  order  to  find  max  J3(/i),  we  must  find  the  h  >  0  that  solves 

h 


HhL\\A\\)  +  1 
rin(e/C) 


-  2k77i  =  0 


such  that 


1 


2kT  <  0. 


(48) 


(49) 


Thln(e/C) 

We  must  have  T  and  h  positive,  and  for  a  numerical  solution  to  be  meaning¬ 
ful,  certainly  e  must  be  very  small  and  in  general  much  less  than  C .  Thus, 
In  {s/C)  <  0.  We  then  require  k  >  0  to  ensure  that  the  second  derivative  of  J3 
is  negative,  which  guarantees  that  the  solution  to  (48)  is  indeed  a  maximum. 


Now  let 


£Pii 


K=-A22T^4t C)’  (50) 

where  A  is  a  free  parameter  that  weights  the  optimization  toward  efficiency 
in  time  or  toward  global  error.  A  better  understanding  of  how  A  affects  the 
variable  step-size  selection  process  can  best  be  explained  by  studying  Table 
4.1.2  and  Table  4.2.2.  By  substituting  this  k  into  (48),  we  find  that  we  must 
solve 

ln(AL||A||)  +  1  +  X2hL\\A\\  =  0  (51) 

for  h  given  an  arbitrary  value  for  A.  In  practice,  we  actually  make  the  substi¬ 
tution  x  =  hL||A||  and  solve 


In  a:  +  1  +  X2x  =  0 


(52) 


X 


for  x.  We  then  compute  h  =  — - ,, .  The  solution  to  this  equation  exists  and  is 

LU II 

unique.  This  is  because  the  h  that  solves  (51)  is  the  unique  global  maximum 
of  the  function  J3,  which  exists  because  of  the  concavity  of  ./3 .  Furthermore, 
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(51)  must  be  solved  numerically  for  A  7^  0;  for  example,  Newton’s  method  or 


Picard  iteration  may  be  used.  For  A  =  0,  notice  that  the  solution  is  h  =  — ,,, 

eL\\A\\ 

which  was  discussed  earlier. 


If  one  is  interested  in  finding  an  equation  similar  to  (51)  for  a  Runge-Kutta 
method  other  than  the  implicit  midpoint  method,  two  things  change  in  (45). 
First,  1 1 A 1 1  will  certainly  change  when  the  method  changes.  It  will  be  necessary 
to  change  the  second  term  as  well.  Suppose  the  method  chosen  has  global  error 
0(Thr).  Then,  (45)  becomes 


Mh) 


hln(hL\\A\\) 

T]n(e/C) 


-  nThr. 


(53) 


Now  we  define 


k  =  —A2 


2T21ii(e/C)  ’ 


and  discover  that  we  must  now  solve 


(54) 


In  a:  +  1  +  X2xr  1  =  0 

for  x  after  again  making  the  substitution  x  =  hL  1 1 A 1 1 . 


(55) 


3.2  Local  Optimization  Using  Picard  Iteration 

If  one  considers  the  analysis  presented  in  the  previous  section,  an  obvious 
question  might  arise.  Is  it  possible  to  repeat  the  process  while  minimizing  local 
error  and  local  computations  instead  of  global  error  and  total  computations? 
It  turns  out  that  the  answer  is  no. 

Instead  of  minimizing  (41),  i.e.  maximizing  (43),  if  we  want  to  minimize  the 
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total  number  of  computations  locally  we  minimize  the  cost  function 

In  {e/C) 


Ji(h)  = 


Equivalently,  we  can  maximize 


Uh)  = 


In  (ftiimil)  ' 


In  (fti||U||) 


(56) 


(57) 


In  (s/C) 

We  also  want  to  consider  the  local  truncation  error.  Just  as  in  the  previous 
section,  we  choose  the  implicit  midpoint  method  for  explanation  purposes. 
The  local  truncation  error  for  this  method  is  0(h3).  Thus,  similar  to  (45),  we 
want  to  maximize  the  cost  function 


In  (fcLPU) 

Mh)  =  Me/Cj  ~ Kh  ■ 


(58) 


First,  we  compute 


—  3  k/T 


dh  h\n(e/C) 

and  set  it  equal  to  zero.  After  doing  so  and  rearranging,  we  get 


(59) 


hs  = 


1 


(60) 


3k  In  {e/C) ' 

Since  e  C  in  general,  (e.g.  in  our  implementations  we  choose  £  =  10-10),  it 
will  certainly  be  the  case  that  In  (s/C)  <  0.  Thus,  (60)  would  imply  that  it  is 
necessary  to  have  k  <  0,  since  it  must  be  that  h  >  0. 


Next,  we  compute 


d2J, 


-1 


dh2  h 2  In  (e/C) 


—  6k/i. 


Since  we  are  trying  to  maximize  (58),  we  need 


d2J3 

dh2 


(61) 


<  0.  However,  if  we 


consider  each  term  in  (61),  we  can  easily  see  that  this  is  impossible.  Since 


In  (s/C)  <  0,  it  must  be  that 


-1 


h?  In  {s/C) 


>  0.  Since  we  have  determined 
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that  we  must  have  k  <  0,  it  must  also  be  that  —6nh  >  0.  Hence,  the  second 
derivative  would  imply  that  the  function  we  are  optimizing  is  actually  concave 
up  and  has  no  local  (or  absolute)  maximum  for  h  >  0. 

The  analysis  above  shows  that  this  method  of  choosing  variable  step-sizes 
(minimization  of  a  particular  efficiency  function)  fails  if  one  considers  local 
truncation  error  and  local  computations  for  the  efficiency  function.  Although, 
there  certainly  are  ways  of  choosing  h  based  on  local  truncation  errors  as  we 
described  in  the  Introduction,  we  simply  cannot  use  local  error  and  computa¬ 
tions  when  using  the  methods  presented  in  this  paper. 


3.3  Lipschitz  Constant  Unknown 

For  most  initial  value  problems  the  Lipschitz  constant  of  /  is  not  easily  acces¬ 
sible.  In  this  case,  an  approach  that  is  slightly  different  than  that  of  Section 
3.1  is  taken  to  optimize  h.  The  idea  in  this  case  is  to  linearize  the  function  / 
at  each  step  of  integration  by  computing  the  Jacobian  of  /.  We  essentially  find 
an  optimal  h  at  each  step  of  the  integration  using  the  analysis  from  Section 
3.1.  The  method  is  described  in  detail  below: 

(1)  Choose  a  value  for  the  parameter  A.  (A  method  for  choosing  A  will  be 
given  in  Section  4.1.) 

(2)  Solve  equation  (52)  for  x  once. 

(3)  At  t  =  0,  let  L  =  ||.D/||,  where  Df  is  the  Jacobian  matrix  of  /,  and 
compute  h  = 

(4)  Perform  one  step  of  integration  using  the  implicit  midpoint  rule. 

(5)  Recompute  L  using  the  new  values  of  the  state  variables,  and  use  this  L 
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to  find  a  new  optimal  h. 


(6)  Repeat  steps  four  and  five  until  the  integration  reaches  t  =  T. 


3-4  Optimization  Using  Newton’s  Method 


We  mentioned  in  Section  2.1  that  one  might  also  use  Newton’s  method  to  solve 
for  the  stage  values  Y^.  Because  of  this,  the  analysis  for  finding  an  optimal 
value  of  h  can  be  repeated  for  Newton’s  method.  Convergence  properties,  and 
hence  convergence  theorems,  are  more  complicated  for  Newton’s  method  than 
for  fixed  point  iteration.  Because  of  this,  one  might  expect  the  analysis  to  be 
a  bit  more  involved  than  that  of  equations  (34)-(51). 

Before  we  begin  looking  at  the  optimization  process  for  Newton’s  method,  let 
us  first  consider  the  following  Lemma. 

Lemma  3.1  If  R  is  an  invertible  n  x  n  matrix,  then 


R 


-1 


< 


\\R\\ 


n—  1 


|  det  R| 


(62) 


Proof.  Let  A;  (r1R')  denote  the  i-tli  eigenvalue  of  RTR.  Since  RT R  >  0,  we 
have  A*  (^RT Rj  >  0  for  i  —  1, . . . ,  n  and  (RT Rj  >  0.  Now  using  a  well-known 
property  of  matrices  and  the  Rayleigh-Ritz  inequality,  we  get  the  following: 


20 


R-1 


=  Amax((i?TR)_1) 

1 

Amin  (WR) 

yiu\  (rtr) 

~  Amin  (RTR)  •  det  ( RTR ) 

r  /  rr,  \  1  n—  1 

Amax  (rtr) 

~  (det  R)2 
_  ||JR||2(n_1) 

(det  R)2 

iW1^ 

y  det  R  J 


(63) 

(64) 

(65) 

(66) 

(67) 

(68) 


Taking  square  roots  on  both  sides  of  the  above  inequality  yields  the  desired 
result.  □ 


Again,  we  consider  solving  (2)  numerically  on  the  interval  [0,  T]  which  is  par¬ 
titioned  by  the  following  sequence  of  points:  {kh}^=0.  Stoer  and  Bulirsch  [10] 
prove  a  theorem  showing  that  Newton’s  method  is  quadratically  convergent. 
We  will  find  an  optimal  choice  of  h  using  the  results  of  that  theorem.  We  begin 
by  defining 


Pk  ~  - ^ - 


(69) 


where  and  7  are  described  below.  First  we  let  C  be  a  convex  subset  of 

]Rns  and  mention  that  the  theorem  from  [10]  and  the  analysis  below  is  valid 
only  in  a  neighborhood,  Sr  (1^°),  of  the  initial  guess  Y®  such  that  Sr  (Y®)  C  C. 
We  then  let 
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®k  '■ 


Yk  -  n° 


k  =  1, . . . ,  K 


/?,-SUp  II DGi^YW^  k  =  1 

yjy  |det  DG  (4_i,  y) | 

~f  =  hZ\\A\\, 


.,K 


(70) 

(71) 

(72) 


where  L  is  a  constant  that  satisfies 


Df(t,  u)  -  Df(t,v ) 


<  L||m  —  v 


(73) 


for  any  vectors  u  and  v  in  Mn's .  Note  that  this  analysis  holds  whether  the 
Lipschitz  constant  of  Df,  L,  is  given  from  the  beginning  or  if  it  is  approximated 
similarly  to  what  was  done  in  the  previous  section;  we  shall  only  require  that 
L  be  a  known  constant  at  the  7-th  step  of  integration. 


It  should  also  be  noted  that  in  practice,  i.e.  in  the  example  we  show  later,  cq, 
and  f3k  will  not  depend  on  k.  For  implementation  purposes,  we  elect  to  choose 
1)°  to  be  the  same  vector  for  all  k.  Hence,  oq,  will  be  the  same  for  each  k.  In 
general,  this  is  not  necessary;  that  is  why  we  define  a  below  for  convenience  of 
analysis.  Actually,  /3k  must  satisfy  DG (4_i,F)_1  <  f3k  for  all  Y  e  C.  We 
then  apply  Lemma  3.1  to  arrive  at  the  definition  given  by  equation  (71).  Now, 
(3k  depends  on  the  step-size  through  the  Jacobian  matrix  of  G ,  i.e.  through 
DG.  Since  this  is  a  variable  step-size  method,  this  implies  that  (3k  should 
be  computed  at  each  step  of  the  integration  using  the  current  step-size.  We 
quickly  found  that  computing  /3k  at  each  step  of  the  integration  makes  the 
process  computationally  inefficient.  Instead,  we  approximate  (3  =  max  /3k 

k 

before  solving  a  particular  example  by  first  solving  the  system  on  a  much 
smaller  time  interval.  As  we  solve  the  problem,  we  keep  track  of  the  current 
values  of  (3k  and  keep  the  largest  value  to  use  as  a  global  constant  to  solve  the 
entire  example.  Putting  all  this  together  and  using  the  theorem  from  Stoer 
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and  Bulirsch,  we  know  that  for  k  fixed,  the  iterations  of  Newton’s  method 
must  satisfy 


F+1  -  Yj 


/  2^-1  ^  2J  —1 
<  akPk  <  apk  , 


(74) 


where  a  =  max  a*,. 

k 


Furthermore,  we  must  assume  (as  do  Stoer  and  Bulirsch)  that  Pk  <  1.  Then 
for  each  k,  it  is  sufficient  to  stop  iterating  Newton’s  method  at  the  smallest 
natural  number  greater  than  or  equal  to  M,  where  M  satisfies  ap2^~l  =  e. 
Solving  for  M,  we  get 


M 


logs 


/In {fikSo^Y 

V  lnP  , 


(75) 


Again,  we  choose  A/,  =  \M~\ .  Also,  we  are  showing  the  analysis  for  the  im¬ 
plicit  midpoint  method  which  has  global  error  0{Th2).  Thus,  we  minimize  the 
following  cost  function: 


Mp) 


Ta/3L\\M  ,  (Hpkea-y\  AkTp^_ 

2 Pk  °g2\  In p  )  +  {a(3L\\A\\Y  ■ 


(76) 


We  now  define 


K  = 


?WL\m3 

''  16  In  2 


(77) 


Then,  after  using  a  little  calculus  similar  to  equations  (45)- (49),  we  find  the 
optimal  h  by  first  finding  the  pk  that  solves 


(in (pksa  x))  1  -  (In pk)  1  -  (ln2)log2  ^ Pk  =  0,  (78) 


and  then  computing 


h  —  2 pk 

k+l  ~  aflL\\A\\ ' 

Because  of  the  complexity  of  (78),  we  must  solve  for  pk  numerically. 


(79) 
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4  Examples 


In  this  section  we  explore  this  variable  step-size  selection  method  for  two 
problems,  the  Lotka-Volterra  model  and  the  Kepler  problem. 


4-1  The  Lotka-Volterra  Model 


For  this  example  we  consider  the  Lotka-Volterra  model  of  a  simple  predator- 
prey  system  from  mathematical  biology  This  particular  example  is  taken  from 
Hairer,  Lubich,  and  Wanner  [1],  Consider  the  following  system: 


it 

u(v  —  2) 

V 

u(l  —  u ) 

/(«,«);  te[0,50]. 


(80) 


In  (80),  u(t)  denotes  the  number  of  predators  present  at  time  t,  v(t)  represents 
the  number  of  prey  present  at  time  t,  and  the  constants  one  and  two  have  been 
chosen  arbitrarily.  This  system  was  integrated  numerically  using  the  implicit 
midpoint  rule.  Since  the  system  is  non-linear  and  the  Lipschitz  constant  of  the 
system  as  a  whole  is  unknown,  we  will  use  the  method  described  in  Section 
3.3. 


This  procedure  was  compared  to  a  fixed  step-size  integration  method  with 
random  step-sizes  chosen.  Two  measures  were  chosen  for  comparison.  The 
first  measure,  T,  was  total  cpu  time  (in  seconds)  for  1000  runs  with  random 
initial  data  uniformly  distributed  on  [0.1, 10].  The  second  measure,  E,  was  the 
maximum  absolute  variation  of  the  numerical  method  from 


I(u,v)  =  In  u  —  u  +  2  In  v  —  v, 


(81) 
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an  invariant  quantity  for  this  system.  The  initial  data  for  the  system  in  this 
case  was  chosen  to  be  [w(0)  u(0)]T  =  [2  6]T. 

We  found  that  for  simple  systems  such  as  (80),  the  numerical  computational 
overhead  in  computing  the  step-size  in  the  optimal  h  method  renders  the 
method  less  useful  than  a  simple  fixed  step-size  method.  After  trying  various 
fixed  step-sizes,  it  was  determined  that  for  1000  runs  with  random  initial 
data,  h  =  0.125  was  the  largest  fixed  step-size  attempted  that  permitted 
convergence.  For  h  =  0.125,  T  =  118.258  and  E  =  0.064.  For  the  optimal  h 
method,  various  values  for  A  were  tried  until  a  comparable  value  for  E  was 
found.  For  instance,  for  A  =  2  we  get  E  =  0.143;  for  A  =  3  we  get  E  =  0.068; 
and  for  A  =  4  we  get  E  =  0.037.  Since  A  =  3  yielded  a  comparable  value  of 
E,  A  =  3  was  chosen  for  1000  runs  with  random  initial  data  and  it  was  found 
that  T  =  195.570. 


Different  results  arise  when  we  try  more  challenging  problems.  Consider  this 
variation  to  the  Lotka-Volterra  problem: 
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(82) 


This  system  has  has  the  same  invariant  as  (80) ,  but  is  very  sensitive  to  random 
initial  data.  For  this  reason  the  initial  data  is  fixed  at  [u(0)  u(0)]T  =  [2  3]r 
for  the  computation  of  both  T  and  E. 


Two  methods  were  chosen  to  solve  for  the  stage  value  y\  which  is  defined 
implicitly  by  (10).  The  first  method  is  the  algorithm  of  Section  2,  which  is 
simply  a  Picard  iteration.  Secondly,  we  used  Newton’s  method  to  compute 
the  stage  value  y\.  The  results  from  this  example  are  given  in  Tables  4.1.1 
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through  4.1.4  and  Figures  4.1.2  and  4.1.3. 


To  compare  the  fixed  step-size  method  to  the  variable  step-size  method,  we 
must  locate  times  that  are  comparable  from  the  tables  and  then  compare 
the  equivalent  error.  For  example,  we  first  notice  that  for  the  fixed  step-size 
h  =  0.1  in  Table  4.1.1,  the  method  took  160.701  seconds  to  solve  the  problem 
using  a  Picard  iteration  to  solve  for  iji.  The  error  involved  for  this  step-size 
was  0.094.  Now  we  look  in  Table  4.1.2  and  find  that  when  A  =  2,  the  problem 
was  solved  in  168.412  seconds,  which  is  about  eight  seconds  longer  than  for 
h  =  0.1.  However,  we  notice  that  the  error  has  been  reduced  to  0.004,  which  is 
about  4%  of  the  error  from  when  h  =  0.1.  We  can  locate  other  instances  similar 
to  this  from  the  two  tables.  For  the  fixed  step-size  h  =  0.08,  the  problem  is 
solved  in  234.427  seconds  using  Newton’s  method  to  find  i]\ ,  yielding  an  error 
of  0.069.  We  compare  this  to  A  =  2  which  was  solved  in  229.851  seconds  with 
an  error  of  0.004.  In  addition,  for  the  fixed  step-size  h  =  0.125  using  Newton’s 
method  to  solve  for  yi,  the  problem  is  solved  in  155.564  seconds  with  an  error 
of  0.084.  We  compare  this  to  A  =  1  in  which  the  problem  is  solved  in  151.649 
seconds  with  an  error  of  0.025. 

In  Table  4.1.2,  when  we  computed  the  stage  value  ij\  using  Newton’s  method, 
we  actually  used  the  variable  step-sizes  determined  by  the  solution  of  (52)  not 
(78).  For  comparison,  we  recreated  these  table  on  another  machine;  this  time 
we  use  (78)  to  compute  the  variable  step-sizes  when  we  are  solving  for  the 
stage  value  y\  using  Newton’s  method.  The  results  are  given  in  Table  4.1.3 
and  Table  4.1.4.  There  are  a  couple  of  interesting  things  to  point  out.  First,  we 
notice  that  for  fixed  step-sizes,  Newton’s  method  works  faster  than  the  fixed- 
point  iteration.  We  must  note  that  the  MATLAB®  code  that  generated  the 
data  for  Tables  4.1.1  through  4.1.4  was  exactly  the  same;  the  only  difference  is 
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that  Tables  4.1.3  and  4.1.4  contain  data  from  a  different  version  of  MATLAB® 
and  ran  on  a  different  machine.  The  main  goal  of  this  paper  is  not  to  compare 
Newton’s  method  versus  Picard  iteration  as  much  as  it  is  to  compare  fixed 
step- sizes  to  variable  step-sizes.  In  this  regard,  the  results  shown  in  Tables 
4.1.3  and  4.1.4  agree  with  those  given  in  Tables  4.1.1  and  4.1.2.  For  example, 
when  we  use  the  fixed  step-size  h  =  0.05,  the  problem  is  solved  in  67.701 
seconds  with  an  error  of  0.031  when  solving  for  the  stage  value  using  a  Picard 
iteration.  Using  the  variable  step-size  method  with  A p  =  3.6,  we  see  that  the 
problem  is  solved  in  67.413  seconds  with  an  error  of  0.003,  which  is  about  90% 
smaller  than  the  fixed-step  size  error  of  0.031.  A  similar  comparison  can  be 
made  with  h  =  0.125  and  A p  =  2. 

The  second  interesting  thing  we  notice  is  evidenced  by  Table  4.1.4  when  we 
use  (78)  to  compute  the  variable  step-sizes  and  solve  for  the  stage  value  iji 
using  Newton’s  method.  We  notice  that  the  problem  takes  longer  to  solve, 
but  at  the  same  time  the  error  is  much  smaller.  We  cannot  compare  the  fixed 
step-size  data  with  the  variable  step-size  data  as  above  because  of  the  large 
difference  in  times.  However  we  can  note  that  when  the  problem  is  solved  with 
A p  =  0.4  the  problem  is  solved  in  27.097  seconds  with  an  error  of  0.087;  when 
the  problem  is  solved  with  XN  =  10,  it  takes  115.689  seconds  to  get  the  error 
down  to  0.001130.  We  can  quickly  compute  that  the  problem  is  solved  77% 
faster  using  Xp,  but  the  error  is  98%  lower  using  A  at- 

The  reason  it  takes  so  long  to  solve  the  problem  using  this  method  is  because 
of  the  variable  step-sizes  determined  by  solving  (78).  Consider  the  plot  shown 
in  Figure  4.1.4.  This  graph  shows  how  far  the  integration  of  the  problem  has 
progressed  at  each  iteration  for  various  values  of  Xp  and  A  at-  Essentially,  the 
graph  also  shows  how  large  (or  how  small)  of  a  step-size  each  method  chooses 
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as  the  integration  proceeds  from  one  step  to  the  next.  In  addition,  we  can 
see  from  the  plot  that  when  we  use  (78)  to  find  the  step-size  update,  the 
step-sizes  determined  are  much  smaller  then  when  (52)  is  used;  hence,  many 
more  iterations  are  required  to  integrate  the  system  from  t  =  0  to  t  =  50. 
We  would  also  like  to  mention  something  about  the  existence  of  a  solution  to 
(78).  Although,  we  have  not  proven  that  a  solution  does  exist,  we  can  plot  the 
function  for  reasonable  values  of  the  parameters.  For  example,  consider  Figure 
4.1.5.  Here  we  plot  the  function  g(p)  where  g  is  the  function  given  by  the  left 
hand  side  of  (78).  This  plot  uses  the  following  values  for  the  parameters  given 
in  (78):  A  =  10,  £  =  10-10  and  a  =  [2  3]T  .  The  function  appears  to  be  quite 
smooth  in  this  region  and  clearly  we  see  that  for  these  parameters,  a  solution 
to  (78)  exists. 


As  one  can  see  from  the  example  above,  inherent  with  this  variable  step-size 
selection  method  is  the  choice  of  the  parameter  A.  We  will  use  the  system  given 
by  equation  (82)  to  explain  how  one  should  choose  an  appropriate  value  of  A 
when  integrating  a  system  that  evolves  over  a  long  period  of  time.  Suppose 
we  are  interested  in  integrating  the  system  described  by  (82)  over  the  interval 
[0,  500]  or  [0, 1000].  First,  we  choose  a  much  smaller  value  for  the  final  time 
of  integration;  in  this  example  that  value  is  T  —  50.  We  then  integrate  the 
system  over  the  interval  [0,  50]  with  a  fixed  step-size  and  at  the  same  time 
with  various  values  of  A.  Essentially,  we  analyze  how  A  affects  this  system 
in  particular,  just  as  we  did  in  the  above  example.  After  we  have  integrated 
the  system  over  the  much  smaller  time  interval,  we  choose  the  value  of  A  that 
works  best  for  this  system  to  integrate  the  system  over  the  entire  time  interval. 
This  process  should  be  done  for  any  system  where  the  length  of  the  interval 
over  which  the  integration  must  be  performed  is  quite  large  when  compared 
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to  the  evolution  of  the  dynamics  of  the  system. 


All  computations,  except  those  that  generated  Tables  4.1.3  and  4.1.4,  were 
done  in  MATLAB®  version  6.1.0.450  Release  12.1  running  in  Microsoft  Win¬ 
dows  XP  Professional  version  5.1.2600  with  an  Authentic  AMD  processor  run¬ 
ning  at  1544  Mhz. 


4-2  The  Kepler  Problem 

This  example,  taken  from  Hairer,  Lubich,  and  Wanner  [1],  is  the  well  known 
example  describing  planetary  motion  discovered  by  J.  Kepler.  We  shall  con¬ 
sider  the  following  two-body  problem: 


<?i  = 
<?2  = 


Qi 


<9?  +  <gfn 
Q2 


(9f  +  9l)3/2’ 
with  the  following  initial  conditions: 

<?i  (0)  =  1  —  e,  92(0)  =  0,  9i(0)  —  0,  92(0)  = 


(83) 

(84) 


'1  +  e 
1  -  e 


(85) 


where  0  <  e  <  1.  To  describe  the  motion  of  two  bodies,  one  of  the  bodies 
is  taken  to  be  the  center  of  the  coordinate  system  and  the  position  of  the 
second  at  any  time  t  is  given  by  the  two  coordinates  cp  (t)  and  ry2(t).  Equations 
(83)- (84)  are  equivalent  to  the  following  Hamiltonian  system: 


H{pi,P2,qi,q2)  = 


T=Pi  *  =  1,2 

(p?+pS)  1 


qi  +  ql 


(86) 

(87) 


where  H{p\,  q±,  q-i)  is  the  Hamiltonian  of  the  system. 
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Fig.  4.1.3.  Error  comparison  for  h  and  A  for  the  Lotka-Volterra  model 
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Before  (83)-(84)  can  be  integrated  using  the  implicit  midpoint  rule,  we  convert 
the  equations  to  a  system  of  four  first-order  ordinary  differential  equations. 
Let  z\  =  qi,  z2  =  qi,  z3  =  q2  and  z4  =  q2.  Define  z  =  [z±  z2  z3  Z4 ]r.  Then 
(83)- (84)  are  equivalent  to  the  following  system: 


Zi 

2(2 

ig 

21 

= 

(4+4f/2 

Z3 

z4 

z4 

23 

L  (*l2+2§)3/2  j 

(88) 


The  above  system  of  equations  was  solved  using  the  implicit  midpoint  rule 
for  t  e  [0,50].  Just  as  in  the  previous  example,  both  a  Picard  iteration  and 
Newton’s  method  were  used  to  compute  the  stage  value  y\ .  The  two  measures 
we  chose  for  this  example  are  very  similar  to  the  those  of  the  previous  example. 
The  first  measure  is  T,  the  total  cpu  time  required  to  solve  the  system  1000 
times  with  eccentricity,  e,  that  is  uniformly  distributed  in  the  interval  [0.4, 0.8]. 
The  second  measure  was  E,  the  maximum  absolute  error  that  the  integration 
deviates  from  the  exact  solution  with  eccentricity  e  =  0.6.  For  this  measure, 
we  considered  the  absolute  error  at  every  step  of  the  integration. 


The  computations  were  performed  on  the  same  machine  as  the  previous  ex¬ 
ample,  and  the  results  are  summarized  in  Tables  4.2.1  and  4.2.2  and  also  in 
Figures  4.2.1  and  4.2.2.  We  compare  the  performance  of  the  variable  step-size 
method  to  the  fixed-step  size  method  exactly  as  we  did  in  the  previous  exam¬ 
ple.  In  Table  4.2.1  we  begin  with  the  fixed  step-size  h  =  0.01.  The  problem  is 
solved  in  84.842  seconds  using  Picard  iteration  to  find  yi ,  giving  an  error  of 
0.113.  We  compare  this  to  A  =  8  in  Table  4.2.2  where  the  problem  is  solved 
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in  81.388  seconds  with  an  error  of  0.067.  Also,  when  the  fixed  step-size  is 
h  =  0.05,  the  problem  is  solved  in  21.190  seconds  using  Picard  iteration  to 
find  y\  and  is  solved  in  29.602  using  Newton’s  method  to  find  y\ .  The  error  in 
both  cases  is  1.581.  We  compare  these  to  when  A  =  2  in  Table  4.2.2.  Respec¬ 
tively,  the  times  are  21.572  seconds  and  29.943  seconds.  The  error  for  these 
two  times  is  0.643. 


4-3  The  Kepler  Problem  Using  Runge-Kutta-Ny strom 


Once  again,  we  consider  the  Kepler  problem  described  by  equations  (83)-(85). 
We  now  solve  the  problem  using  a  Runge-Kutta-Nystrbm  (RKN)  method. 
The  method  we  use  is  the  two-stage  Gauss  method  from  [3]  which  is  the  RKN 
method  that  is  induced  from  the  fourth  order  Gauss  method  (a  Runge-Kutta 
routine)  found  in  [1] .  For  notation,  we  will  refer  to  this  integration  routine  as 
Gauss4.  The  constants  from  equations  (28)- (30)  are 
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The  analysis  leading  to  equation  (51)  used  the  Lipschitz  constant  from  the 
Runge-Kutta  methods,  /iL||A||.  Since  we  are  actually  implementing  a  RKN 
method  to  solve  this  system  now,  we  have  to  repeat  that  analysis  process  with 
the  Lipschitz  constant  from  the  RKN  methods,  h2L||A||.  Following  exactly  the 
same  process  given  in  equations  (34)-(52)  and  using  the  substitution 


x  =  h2L\\A\\ , 


(91) 
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we  find  that  to  find  the  optimal  h  we  must  solve 


2  +  In  x  +  X2x 3/2  =  0 


(92) 


for  x  and  then  solve  for  h  using  (91). 

Tan  [5]  solves  the  Kepler  problem  using  Gauss4  as  a  Runge-Kutta  routine  (not 
in  its  induced  RKN  form)  with  a  fixed  step-size  of  ^  for  500,000  iterations. 
Just  for  comparison  we  did  the  following: 

(1)  We  solve  the  problem  using  the  fixed  step-size  ^  on  the  interval  [0,50] 
and  determine  the  maximum  absolute  variance  from  the  exact  solution. 

(2)  We  compute  the  total  cpu  time  taken  to  solve  the  system  using  this  step- 
size  for  500,000  iterations  as  Tan  does. 

(3)  We  solve  the  system  on  the  interval  [0,50]  using  the  optimal  step-size 
selection  method  to  determine  the  value  of  A  that  gives  us  a  comparable 
error  as  determined  in  step  one. 

(4)  We  solve  the  system  on  the  entire  interval  using  this  value  of  A  and  our 
variable  step- size  selection  method. 

When  we  used  h  =  ^  to  solve  the  system  on  the  interval  [0,50],  we  found 
the  maximum  absolute  variance  from  the  exact  solution  to  be  0.0045.  Then 
using  this  fixed  step-size,  we  found  that  it  took  101.577  cpu  seconds  to  solve 
the  problem  taking  500,000  iterations.  We  then  determined  that  a  comparable 
error  on  the  interval  [0,50]  was  achieved  for  A  =  55;  that  error  was  0.0046. 
When  we  solve  the  problem  using  our  variable  step-size  selection  method  and 
choose  A  =  55  to  a  final  time  of  T  =  50°o°Q7r ;  we  find  that  the  system  is  solved 
in  61.461  cpu  seconds,  which  is  39.5%  faster  than  using  the  fixed  step-size  of 
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In  addition  to  the  above  comparison,  we  also  compared  the  fixed  step-size 
method  to  the  variable  step-size  method  for  various  values  of  h  and  A  for  this 
example.  The  results  of  this  comparison  are  given  in  Tables  4.3.1  and  4.3.2. 
We  only  used  Picard  iteration  to  solve  for  the  stage  values.  We  used  exactly 
the  same  two  measures  as  we  did  in  the  previous  section,  when  we  solved  the 
Kepler  problem  using  the  implicit  midpoint  rule.  When  we  solve  the  system 
with  the  fixed  step-size  h  =  0.05,  the  solution  is  found  in  198.247  seconds  with 
an  error  of  0.0050.  Comparatively,  we  found  that  for  A  =  100,  the  system  is 
solved  in  193.291  seconds  with  an  error  of  0.0036,  which  is  about  28%  lower 
error.  We  also  found  that  when  the  system  is  solved  with  the  fixed  step-size 
of  h  =  0.06,  it  took  179.498  seconds  with  an  error  of  0.0071.  We  compare  this 
to  the  variable  step-size  with  A  =  90,  which  took  only  182.738  seconds  to  once 
again  get  an  error  of  0.0036,  which  is  about  49%  lower  error.  In  addition  to 
these  positive  results,  we  did  not  notice  a  large  decrease  in  the  error  using 
the  variable  step-size  method  as  we  increased  A.  When  we  used  the  implicit 
midpoint  rule  in  the  previous  two  examples,  we  noticed  a  steady  decline  in  the 
error  as  A  increased.  In  this  case,  we  actually  noticed  that  the  error  only  goes 
from  0.0036  when  A  =  90  to  0.0032  when  A  =  1000.  We  point  out  that  when 
A  =  1000,  the  system  is  solved  in  625.772  seconds  which  is  about  33%  faster 
than  the  927.805  seconds  it  took  with  the  fixed  step-size  h  =  0.01,  where  the 
error  was  only  0.0029. 


All  computations  were  done  in  MATLAB®  version  6.5.0.180913a  Release  13 
running  in  Microsoft  Windows  XP  Professional  version  5.1.2600  with  an  x86 
Family  15  Model  2  Stepping  7  Genuinelntel  processor  running  at  2392  Mhz. 
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5  Conclusion 


The  collection  of  implicit  numerical  integration  routines  is  vast  to  say  the  least. 
Often  times  one  routine  is  chosen  over  another  to  improve  either  efficiency  or 
accuracy.  In  this  paper,  we  have  shown  that  it  is  possible  to  wisely  choose  a 
variable  step- size  for  these  integration  schemes. 

For  linear  ordinary  differential  equations  or  equations  in  which  the  Lipschitz 
constant  for  the  function  /  is  known,  the  task  becomes  quite  simple  as  the 
optimal  value  of  the  step-size  will  not  change  from  one  step  of  the  integration 
to  the  next.  But,  if  we  are  dealing  with  more  complicated  non-linear  differential 
equations,  we  can  still  choose  an  optimal  time  step  at  each  step  of  integration  of 
the  system.  As  we  have  shown,  this  process  often  involves  solving  a  non-linear 
equation  numerically.  Because  of  this,  the  computational  overhead  in  using 
this  optimal  step-size  routine  seems  to  be  too  much  for  solving  differential 
equations  in  which  the  function  /  is  quite  simple.  However,  our  results  have 
shown  that  this  is  consistently  not  the  case  when  /  is  a  complicated  function 
as  we  describe  below. 

Tables  4.1.1,  4.1.2,  4.1.3,  4.1.4,  4.2.1  and  4.2.2  clearly  show  that,  for  compa¬ 
rable  integration  times,  the  variable  step- size  selection  method  presented  in 
this  paper  drastically  reduces  the  global  error  in  the  solution  of  the  problem. 
For  the  Kepler  problem,  we  found  that  for  comparable  execution  times,  the 
error  was  reduced  41%  to  59%  when  the  variable  step-size  method  is  used.  In 
the  Lotka-Volterra  example,  we  found  that  for  comparable  execution  times, 
the  problem  is  solved  with  the  error  reduced  70%  to  96%  when  we  use  the 
variable  step-size  method.  From  studying  the  tables  we  may  choose  A  so  that 
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execution  times  are  comparable,  in  which  case  the  variable  step-size  method 
noticeably  reduces  error  as  evidenced  from  the  above  discussion.  However,  A 
may  also  be  adjusted  to  find  comparable  errors  between  the  fixed  step-size 
and  variable  step-size  methods.  When  you  do  this,  one  notices  that  the  time 
required  to  achieve  a  comparable  error  for  the  fixed  step-size  is  much  larger. 

We  must  point  out  that  this  optimal  step-size  selection  process  is  dependent 
upon  the  scheme  being  used  and  we  have  concentrated  on  Runge-Kutta  and 
Runge-Kutta-Nystrom  methods.  It  should  not  be  too  difficult  of  a  task  to 
adapt  this  process  to  the  ever  growing  collection  of  implicit  integration  rou¬ 
tines. 
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Table  4.1.1 


Fixed  step-size  (Lotka-Volterra  model) 


h  -*■ 

0.05 

0.08 

0.1 

0.125 

0.2 

0.25 

T  (Picard) 

239.995 

181.081 

160.701 

159.590 

DNC 

DNC 

E  (Picard) 

0.031 

0.069 

0.094 

0.084 

DNC 

DNC 

T  (Newton) 

354.380 

234.427 

187.720 

155.564 

104.811 

93.294 

E  (Newton) 

0.031 

0.069 

0.094 

0.084 

0.228 

0.228 

Table  4.1.2 

Variable  step-size  (Lotka-Volterra  model) 


A  -► 

0.25 

0.4 

0.5 

0.75 

1 

2 

T  (Picard) 

113.834 

119.692 

118.531 

124.759 

127.323 

168.412 

E  (Picard) 

0.115 

0.087 

0.078 

0.048 

0.025 

0.004 

T  (Newton) 

117.586 

124.294 

127.428 

140.652 

151.649 

229.851 

E  (Newton) 

0.115 

0.087 

0.078 

0.048 

0.025 

0.004 
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Table  4.1.3 


Fixed  step-size  (Lotka-Volterra  model) 


h  — > 

0.05 

0.065 

0.08 

0.1 

0.125 

T  (Picard) 

67.701 

57.018 

50.742 

46.178 

45.596 

E  (Picard) 

0.031 

0.050 

0.069 

0.094 

0.084 

T  (Newton) 

67.801 

53.260 

43.990 

35.574 

29.636 

E  (Newton) 

0.031 

0.050 

0.069 

0.094 

0.084 

Table  4.1.4 

Variable  step-size  (Lotka-Volterra  model) 


Xp  —> 

0.4 

0.6 

0.8 

1 

2 

3.6 

T  (Picard) 

27.097 

28.377 

28.761 

33.902 

43.595 

67.413 

E  (Picard) 

0.087 

0.066 

0.044 

0.025 

0.004 

0.003 

Xn  — ► 

8 

10 

12 

14 

16 

20 

T  (Newton) 

192.345 

115.689 

138.058 

134.097 

140.481 

149.867 

E  (Newton) 

0.000459 

0.001130 

0.000841 

0.000891 

0.000819 

0.000734 
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Fig.  4.1.4.  Comparison  of  the  variable  step-sizes  determined  by  solving  (52)  versus 
(78) 


Table  4.2.1 

Fixed  step-size  (Kepler  problem) 


h  -*• 

0.01 

0.05 

0.1 

0.125 

0.2 

0.25 

T  (Picard) 

84.842 

21.190 

13.299 

11.858 

DNC 

DNC 

E  (Picard) 

0.113 

1.581 

2.038 

2.483 

DNC 

DNC 

T  (Newton) 

137.506 

29.602 

16.503 

14.391 

DNC 

DNC 

E  (Newton) 

0.113 

1.581 

2.038 

2.483 

DNC 

DNC 
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Fig.  4.1.5.  Plot  of  p  versus  g(p),  where  g  is  given  by  the  left  hand  side  of  (78) 


Table  4.2.2 

Variable  step-size  (Kepler  problem) 


A  -*• 

i 

2 

4 

6 

8 

10 

T  (Picard) 

20.199 

21.572 

37.234 

58.505 

81.388 

111.440 

E  (Picard) 

1.250 

0.643 

0.295 

0.133 

0.067 

0.037 

T  (Newton) 

22.062 

29.943 

50.602 

79.263 

114.673 

157.194 

E  (Newton) 

1.250 

0.643 

0.295 

0.133 

0.067 

0.037 
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Fig.  4.2.2.  Error  comparison  for  h  and  A  for  the  Kepler  problem 
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Table  4.3.1 


Fixed  step-size  (Kepler  problem  using  RKN) 


h  -»• 

0.01 

0.025 

0.05 

0.06 

0.1 

T  (Picard) 

927.805 

341.167 

198.247 

179.498 

DNC 

E  (Picard) 

0.0029 

0.0030 

0.0050 

0.0071 

0.0508 

Table  4.3.2 

Variable  step-size  (Kepler  problem  using  RKN) 


A  -*• 

90 

100 

250 

500 

1000 

T  (Picard) 

182.738 

193.291 

248.634 

352.381 

625.772 

E  (Picard) 

0.0036 

0.0036 

0.0032 

0.0032 

0.0032 
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